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ABSTRACT 


Three  dimensional  off  design  flow  fields  are  calculated  for  stream 
Mach  numbers  in  the  range  1.3  to  4.0  and  corresponding  to  attached  and 
detached  shocks  at  the  leading  edges  of  a reentrant  pyramidal  waverider 
geometry.  The  MacCormack  shock  capturing  version  of  the  Lax-Wendroff 
finite  difference  technique  is  used  with  grids  chosen  to  align  with 
surface,  symmetry,  and  approximate  shock  traces  in  the  transverse  plane. 
Separate  natural  grid  systems  are  defined  for  the  compression  and  ex- 
pansion regions,  and  an  alternating  region  algorithm  is  used  in  combination 
with  a sequential  transfer  of  the  edge  region  boundary  conditions.  The 
latter  are  derived  from  overlapping  portions  of  the  computational  grids 
as  integration  proceeds  axially  to  an  asymptotic  conical  field.  Equi- 
valent attached  shock  cases  result  from  either  of  two  approaches:  the 
alternating  region  algorithm,  or  a consideration  of  solely  the  compression 
region  with  uniform  unknown  conditions  assumed  near  the  edges.  For 
detached  shock  cases  overall  lift  and  drag  coefficients  exhibit  smooth 
variations  between  the  attached  edge  and  detached  apex  limits. 
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1.  INTRODUCTION 


Waverider  configurations  have  been  of  interest  for  some  time  for 
supersonic  flight  since  the  concept  offers  an  implied  control  over  con- 
fined pressure  fields  for  practical  three  dimensional  bodies.  Concep- 
tually an  inverse  approach  is  adopted  to  determine  geometric  surfaces 
consistent  with  fields  derived  from  specified  shock  surfaces.  Typically, 
but  not  necessarily,  the  shock  surface  is  of  a simple  kind  (planar  or 
circular  conic,  e.g.),  and  the  associated  body  corresponds  to  that 
specific  (on-design)  flow  field.  Nonweiler's  original  suggestion  ^ 
was  a delta  planform,  caret  cross  section  wing  of  finite  thickness  and 

anhedral  (Fig.  1)  which  matched  a planar  shock  surface  extending  between 

(2)  (3) 

the  leading  edges  on  the  compression  side.  Gonor'  ' and  Maikapar  1 
considered  an  integral  number  of  similar  components  arranged  circum- 
ferentially so  as  to  form  a star  section,  right  reentrant  pyramid. 

A number  of  studies  have  indicated  that  generalized  waveriders  under 

on-design  conditions  offer  advantageous  lift  to  drag  ratios  and  drag 

(4  5) 

reductions  for  equivalent  volumes.'  ’ ' 

The  essential  geometric  features  of  the  waverider  cross  section 
are  the  external  and  internal  corners  at  the  side  edge  and  midspan 
(or  reentrant  rib)  locations.  The  same  features  are  common,  of  course, 
to  edges  and  junctions  of  many  geometric  components  of  aeronautical 
interest  and  in  particular  for  high  speed  engine  inlets  of  a truncated 
pyramidal  kind. 
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While  on-design  field  evaluations  are  relatively  simple,  the 
off-design  condition  must  be  anticipated  in  most  instances,  and  is  in 
fact  assured  if  viscous  influences  were  to  be  taken  into  account.  The 
shock  surfaces  are  then  unknown,  rotational  fields  are  certain  to  result, 
and  shock  detachment  may  occur.  Several  studies  have  explored  the 
off-design  field  on  the  basis  of  integral  methods,  linear  departures,  and 
numerical  methods,  all  for  attached  shock  cases. The  detached  leading 
edge  case  was  considered  on  a corrected  Newtonian  theory  basis  and 
free  flight  data  for  modified  (blunt  edge)  caret  sections  is  available. 

The  present  objective  is  the  application  of  a shock  capturing  finite 
difference  method  to  caret  waveriders  for  such  off  design  conditions 

i ' ■ 

I 

including  detachment  from  the  leading  (side)  edges.  A second  order, 
explicit  algorithm  is  employed  with  grid  systems  defined  explicitly  for 
corner  section  configurations  and  approximate  alignment  with  average 
shock  trace  locii.  For  the  detached  shock  cases  an  alternating  compression/ 
expansion  region  algorithm  is  used  during  the  convergence  process  following 
the  "time  like"  axial  coordinate  direction. 
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2.  GEOMETRY,  ON-DESIGN  BASIS,  OFF-DESIGN  LIMITS 

The  caret  wing  geometry  may  be  expressed  in  terms  of  the  span,  2s, 
and  the  on-design  flow  deflection,  6,  and  shock,  Bq,  angles  (Figure  1). 
For  Mach  numbers  other  than  MQ  the  section  angle  normal  to  the  leading 
edge  controls  shock  attachment  if  the  shock  surface  remains  attached  at 

i ■ i 

the  apex.  The  normal  section  angle  is 


JN  = cos 


-1 


cosBd  cos6  si n2ipu  - sinBp  sinS  cos2^u 


sini(ju[cos2BD  - dos2<jju  cos2(3D-6)] 


1/2 


; , in!-.' 

where  the  upper  surface  apex  angle  is  defined  by 


*u  = cos"1  I 


COSBr 


[1  + (s*coseD)2]1/2 


J 


(2.1) 


(2.2) 


and  similarly 


' Po  = COS-1  I 


cos^u  cos6(l  + tanBD  tan 6) 


Here  s*  is  the  local  semispan  in  units  of  axial  distance  from  the 
nose,  or  equivalently 


tan  3r 


t*  = 


tan  3, 


B 


(2.3) 


(2.4) 


in  terms  of  the  on-design  shock  angles  in  sections  parallel  and  trans- 
verse to  the  mainstream. 
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The  on-design  condition  is  that  for  which  identical  disturbance 
fields  result  for  oblique  shocks  based  on  either  (Mg,  6)  or  (M^,  6^)  pairs. 

The  normal  section  angle,  6N,  may  be  greater  or  less  than  6 depending 
upon  the  relative  span  (Figure  2),  and  for  attached  shocks  the  normal 
component  of  Mach  number  at  the  leading  edge  is 

i 1 i 

MN  = Moo  sin^u  = Ma>  Sln  [tan_1(s*2  + tan2^)^2]  (2.5) 

Thus  the  geometry  defines  a wide  class  of  swept  wing  flow  field 
behavior  in  the  region  "adjacent"  to  the  leading  edge  if  subject  to 
attached  oblique  shock  conditions.  Typical  parametric  variations  for 
geometries  based  on  Mg  = 2 and  4 are  listed  in  Table  1.  Figure  3 
illustrates  the  span  influence  on  the  presence  of  weak  or  strong  shocks 
at  the  side  edge  for  the  Table  1 cases  applied  to  selected  Mg  and  Bg 
levels. 

A useful  bound  to  the  disturbance  region  when  shocks  are  attached 
is  provided  by  reinterpreting  the  geometry  of  eq.  (2.1).  With  6^  and 
6 replaced  by  and  gQ,  i.e.  by  the  shock  angle  normal  to  the  leading 
edge  and  the  angle  (Figure  1)  formed  by  the  upper  ridge  line  and  the 
intersection  of  the  shock  surface  with  the  symmetry  plane, 

cos gn  cose  sin2^,,  - singn  singn  cos2^,, 
cos  gN  = — 5 5 y D O — U (2.6) 

Simj>u  [cos 2(3q  - cos2^u  cos2(eD-60)]  ' 

This  is  an  implicit  relation  for  the  farthest  possible  extint  of  that 
shock  surface;  i.e. 

60  • so(6r  bdiS*>  (2-7) 

1 
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The  on-design  value  of  $N  follows  from  eq.  (2.6)  on  taking  3Q  = 8D; 
there  results 


cos(3N) 

N D 


r 


s* 

+ tan26D 


cos3d 


i > i 

Typical  values  for  3^  and  8Q  with  departure  of  M^  from  the  design 
value  are  presented  in  Tables  2 and  3. 

The  design  condition,  geometry  and  stream  Mach  number,  also 
can  be  viewed  as  the  configuration  for  which  the  disturbance  is  im- 
parted no  lateral  momentum.  The  flow  downstream  of  the  shock  surface 

(il!<  i 

proceeds  parallel  to  the  geometric  symmetry  plane.  A measure  of 
off-design  is  therefore  the  turning  relative  to  that  plane,  since  the 
symmetry  condition  then  implies  that  some  corresponding  compression 
or  expansion  must  appear  in  order  to  redirect  the  flow.  The  flow 
direction,  t,  from  the  leading  edge  and  in  the  plane  of  the  lower 
(compression) surface  follows  from  the  component  Mach  numbers  after 
the  attached  shock  (Figure  4) 

1 + *T~  (MNsinBfp2 

mn2  = 

2 sin2(8N-<SN)[y(MNsin8N)2  - 


(2.8) 


mt2 


(M^cosi^)2 


1 , 2(y-D(MN  sinVM~ 
(Y+l)2  (Mn  sin8N)2 


y(Mn  sinBN)2  + 1] 


(2.9) 


1 

i 


I 


I.e., 


tan  t = 


(mn  sin3N)2+  5 


6(Mn  sineN)  sin(BN-6N)Moocos^ 


On-design  this  reduces  to  td  = and  (td-  t)  > 0 implies  flow 
towards  either  the  side  edge  or  the  symmetry  plane.  Table  3 includes 
some  typical  turning  levels  and  the  corresponding  lateral  component  of 
velocity. 

. , .4  .l.i  I • I 
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3.  DESCRIPTIVE  FIELD  EQUATIONS  AND  BOUNDARY  CONDITIONS 
Conservation  Equations 

The  governing  equations  for  inviscid,  steady  flow,  of  a perfect 
gas,  with  pressure  and  density  in  units  of  undisturbed  stream  stagnation 
values  and  velocity  components  in  units  of  the  maximum  adiabatic  velocity, 
expressed  in  Cartesian  coordinates  are 


Ex  * Fy  + G2  - ° 


where 


pu 

PU 

■pw,  , 

kp  + pu2 

F = 

puv 

. 6 = 

puw 

puv 

9 * 

kp  + pv2 

» u 

pvw 

puw 

pvw 

kp  + pw2 

and  k = (y-l)/2y.  The  adiabatic  energy  equation  then  specifies 

P = P (1  - q2) 


For  attached  nose  shocks  an  axial  scale  length  is  absent  and  a 
conical  coordinate  system  is  appropriate.  If  (Figure  5) 


C 

n 

Z 


= £nx 


XL 

s 


-JL.  ) 


f ( 


XS 

z 


XS 


then  the  (x,y,z)  to  ( c,n,C)  transformation  implies 


(3.1) 


(3.2) 


(3.4) 


(3.5) 
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> 


_3_  = 1 3 . 3 jL  _ I JL 
3x  _ x 9n  x 3£ 


3 1 3 

~Sy  Its*  3n 


_3_  1 3 ! 

3z  " XS*  3? 


and  so  eqs ,(3.1 ) become  the  weak  conservation  form 


E?  + Fn  + + H = 0 

E = E , F = j*  - nE 

H = 2E  , 6 = |*  - CE 


(3.6) 


(3.7) 


(3.8) 


Boundary  Conditions 

Say  the  body  and  shock  surfaces  are  given  by 

B(r,  a,  0)  = 0 , S (r,  o,  0)  =0  (3.9) 

in  a spherical  (r,  o,  0)  coordinate  system.  A conical  body  geometry 
may  then  be  taken  in  the  form 

Bc(a,  0)  = a - f ( 0 ) = 0 (3.10) 

with  a unit  surface  normal  n = (n^.,  n0>  n0)  given  by 
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(H)T  + (IH)j  + ( — — -^)k 

VB  _ v3r;  V 3 o'J  Vsina  39;k 


n = “JvbT 


r(ll)2  + (1H  )2+  ( — I — H)2]1^2 
*-'3r'  'r  3a  ' Vsino  39'  J 


(sina)j  - (f 1 )k 


(sin2a  + f12) 


TfZ 


. 1 i 


The  unit  surface  tangent  t = (tr>  ta>  t0),  orthogonal  to  7 and  n in 
a positive  9 sense,  is 


t = TJLS  = 


i x n 


= ■ V + no"k  = (°*  "ne»' na)  ■ 


The  velocity  vector  q = (u  , v , w ) has  components  at  the  surface 

c c c 

which  are 


v sin  a - w f* 

q.,  = (q  • n)n  = — - ^TFT  " 

N (sin2a  + f,2)1/2 

v f 1 + w sin  a 

qr  = U_7  + ( q * t ) t = u 7 + yTT  t 

^ c c (sin2a  + f'2)  1/2 


and  for  q^  = 0 it  follows  that 


w. 


sin  o 
f* 


(3.11) 


(3.12) 


(3.13) 


(3.14) 
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For  the  case  of  a planar  conical  surface  such  as  in  Figure  1 the 
eq.  (3.10)  description  implies 


f (0 ) = tan-1 


$*  tan6  sing 


i. 


tangjj  sin(e-g^)  + tan6  cos0  sing^  j 


and 


f'(0) 


si  no  cosa 


tan6  sin0  singQ  - tangncos(0  - g.) 

—P -11 P_ 

tan6  cos0  singg  + tangp  sin(0  -g&) 


sin  o cos  o 
tan(0  - gg  + 6^) 


(3.15) 


(3.16) 


Here  (0-gQ  +6a)  is  the  angle  in  the  transverse  (y,z)  plane  between 

P P 

the  surface  and  the  0 = constant  plane.  Thus  eq.  (3-14)  becomes 


and  implies  that 


tan(0  - Bg  + 6p) 
cos  a 


(3.17) 


< 

W 

v,  1 

<e  - VV  * . 

[0  cw 

t r corresponds  to 

cw 

• = 0 when 

w_ 

lie  in  the  surface  plane 

2 cw  . 

► 

cw 

(3.18) 

The  normal  and  tangent  expressions  eqs.  (3.11),  (3.12)  also  are  valid 
for  a conical  shock 


Sc(o,0)  = o$  - f$(0)  = 0 (3.19) 

j 


n 


and  In  the  conical  system  the  components  of  the  free  stream  velocity 


= (q^coso)!  - (q^sino)] 


(3.20) 


• ■ • 


q sin  a 

'00  c 


(sin’0$  * "• 


q f'  si  no 
00  s s 


5 T = (qcosa  )T  - 5 (f‘  J + sino  k) 

^ ^ s sin2o  + f’2  5 S 

s s i 


(3.21) 


The  shock  angle  is 


3(0)  = sin' 


• 1 1 ^°°N  I 


. -i 
sin 


sin2cr 


(sin2o$  + fj2)1/,Z 


(3.22) 


and  the  post  shock  velocity  components  follow  from  the  usual  shock 
relations  and 


qs  = qs  + V 

S SN  I 


qc  = - qc  nc 

SN  SN  S 


(3.23) 


I.e., 


(Y-l)(MoosinB)2  + 2 

q°°N  (Y+1)(M  sing)2 


(3.24) 


I 
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u = q cosa 
s “ s 


tt  sin2o_  + f'2 


Ps 

r1  - 1 


v , i o » u u • » 0 U 

v = -q  sina  — + = v 1 

sin2as  + fj2  (Y+l)(Moosinas)2  , °°s  y(Moosinos)2 


(3.25) 


V 4r  ~Z.  T 


sin2a. 


q“°~  pJ  f,  ] l'/2 


s I (M^sin  a$)2  sin2a$  + f^  2 


Y M 2 P sir,2as 


where 


, , ■ *•  i .i  ! ■ 1 


P = JLt I_  fi  + Yii 

2yM  2 lP-  Y+1 


Boundary  Conditions  for  Cartesian  Conical  System 


(3.26) 


For  the  modified  Cartesian,  conical,  system  of  eqs.  (3.5)  - (3.8): 


Bc  = K - g(n)  = 0 


(3.27) 


s*(ng'-5)T  - g'j  + k 
(ts*(ng'-5)]2  + g12  + l)1 


<V  nn’  "{> 


The  surface  boundary  condition  q^  = (q*n)n  = 0 then  implies 


„ . w 

nt  + U n£ 


s*(ng'  -0  + J 


n Jvt 


(3.28) 


l 


i 


r 
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Also  the  tangent  component  is  the  surface  velocity  vector 

wT)  E \ = u„ 


f n + (-)  nr  1 

i - c „■  u » s 3 + (77)  e 


w 


(3.29) 


Specializing  to  a planar  conical  surface  as  in  Figure  1,  eq.  (3.27)  is 
specifically  (Figure  5) 


i = 


tanBp  - tan6 


. tan6 


(tane  ) n 


[lower  surface] 
'[tipper  surface] 


E.g.,  for  the  lower  surface  eqs.  (3.28),  (3.29)  are: 


(-) 

vu 


0 - tan6 
w 


w 


(qT) 


= u. 


w 


! (*)  - tanJ 

t ♦ — a 


gj 


3 + 0 R 

w 


(3.30) 


(3.31) 


At  a shock  surface  = gs(n)  the  shock  normal  is  of  the  same  form  as 
in  (3.27).  Now,  however. 


%>  = “co1 


(3.32) 
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4.  FINITE  DIFFERENCE  NUMERICAL  SCHEME 


Predictor  Corrector  Forms 

The  solution  of  the  descriptive  system,  eqs.  (3.7)  was  carried  out 

using  a finite-difference  technique  with  the  computation  advanced  along 

» ■ 1 

the  time-like  axial,  coordinate  until  steady  (axially  uniform) 

conical  conditions  were  achieved.  The  specific  predictor/corrector 

algorithm  was  the  alternating  direction,  explicit  difference  scheme  of 

(12) 

second  order  accuracy  as  suggested  by  MacCormack  ' . 

If  the  differential  description  is  (eqs.  (3.7)) 

j-u  !•  i 

= -Fn  - - R 

and  (n,  j,  k)  are  the  indices  for  grid  lines  in  the  (c,  n»  £)  directions 
then  the  predictor  provides  the  intermediate  state  vector 


(4.1) 


£n+l  = fn  _ AC/pn 
j,k  j»k  Arrrj+.l,k 


FL  > 


I (SL  ♦! 
"L  * 


Gi  i.  > 

J » k 


(4.2) 


and  the  corrector  provides 
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and  from  eqs.  (3.2)  and  (3.8) 

F" 

J>k 

■ f<eL> 

5o".k 

■ E'E".k>  > ■ ’ 

(4.5) 

■ n<Ej".k> 

A solution  is  converged  upon  achieving  constant  e!?  . after  advancing 

J > K 

n a sufficient  number  of  times. 

• i • V-  i.i  ! * . ' 

Grid  Systems 

For  planar  surface  waveriders  the  natural  grid  system  is  one  that 
conforms  to  the  surface,  symmetry,  and  shock  traces  in  the  (n9C)  plane. 
The  shock-capturing  capability  of  the  MacCormack  algorithm  and  the 
unknown  location  of  the  jump  distribution  suggests  the  quadrilateral 
(APBC)  bounded  grid  formed  by  new  coordinates  (x»^)  which  are  defined 
as  (Figure  6) 


* 


i-Lii 

nR  - n 


(4.6) 


These  are  essentially  angular  measures  from  the  surface  and  symmetry 
planes  for  the  lower  (compression)  side.  Approximate  grid  alignment 
with  shock  boundaries  then  follows  from  choices  for  the  (nR,£R)  and 


(q  , 5 ) centers  and  (y  , ip  ).  A procedure  will  be  outlined  below. 
' u u ^nax  max 
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•Inverting  eqs.  (4.6) 


f*u  “ 5R  ' ^nR  , 
n = x(~  > 

i - n 


(4.7) 


e = + ^tnR  - x5y)  , 

1 - 1|>X 


and  furnishes  the  actual  field  coordinates. 

The  transformed  description  in  (c,  x»  4>)  space,  where  x(j)  and  are 
the  normalized  parameters 


S-  .1.!  ! ' > 


x - xO) 


x(V  - xO) 


ip  = 


\p  - ii»(l) 


*(im)  “ MD 


for  each  of  the  two  grid  directions,  follows  from  eqs.  (3.7)  on 

i 

application  of 


_3_ 

3C 


_9_ 

3C 


_3_ 

3n 


r 


3 

35 


n(x(jm)-x(D)  3x  (5-5R)Mim)  -*0))  3* 

y2  3 ^ 3 


— + 


n(x(jm)  - x(D)  3x  U - 5R)(t(im)  - MD)  W 


C4.8) 


(4.9) 
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There  results 


E*  + F*  + G*  + H*  = 0 (4.10) 

C X * 


with  the  revised  state  vectors 


E*  = E 

F*  = X(L+_XG1_ 

n(x(Jm)  - xO)) 

G*  = V + g) 

(5  - eR)(ip(lro)  - ip(l)) 

*[x(Jm)  " x(T ) 3 F*  + xlXO  - 4'0)]G*  XG 

H*  = H + 55 s : - TT-= 


(4.11) 


An  analgous  grid  system  for  the  upper  (expansion)  side  of  the 

* **  . « 

geometry  is  shown  in  Figure  7 and  consists  of 


= £ - nta;i8B 

Or,  on  inverting 

x(su  - O 

n = 

1 + XtanBg 

= X^u  tan^  + ipe 
1 + Xtan30 


(4.12) 


(4.13) 
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r 


Again  normalizing  the  grid  level  as  in  eqs.  (4.8),  the  transformation  to 
(C»  x.  5<e ) space  makes  use  of 


_a_ 


_9_ 


tan  8 


B 


n(x(jm)  “ x(D)  3x  ^e(im)  - ^(1)  W, 


re 

» i 


n(x(jm)  - x(D)  3x  *e(iro)  - 4>eO)  a^e 


and  results  in 


with 


; , ■ i-  i.i  ! i I 


E*  + F*  + G*  + Ho  = 0 
C X %> 


(4.14) 


(4.15) 


E*  -E 


Fe*  = 


Ge*  = 


X (F  + xG) 

n(  x (Jm)  - x(D) 

i 

G - tanBD  F 

D 

*e(1m>  ' *e(,) 


(4.16) 


(1+2  xtanBD)6  + tan8R  F 

H*  = R - [ ? ■? 

e E - \b 

su  ^e 


Grid  Parameters 


The  grid  systems  permit  adjustment  of  the  i and/or  j locii  so  that 
some  mesh  line(s)  approximately  parallel  any  expected  shock  discontinuities 
for  Mro _<  Mp.  Near  on-design  a shock  lying  along  S = corresponds  to  assuming 
SR  ~ for  example.  Off-design  attached  shocks  correspond  to  nR  - nT 
with  ip  = constant  matching  either  the  oblique  shock  trace  near  the  edge,  or 
linearly  approximating  the  entire  shock  trace  over  the  entire  0 £ n £ 1 
interval . 

For  detached  shocks  from  the  side  edge,  choices  for  nR  and  can  be 
made  such  that  the  grid  line  approximates  the  normal  to  the  surfaces, 
and  simultaneously  the  ip  = 0 mesh  line  approximates  the  central  (n  « 1) 
shock  location.  From  eq.  (4.6) 


XH  (^-tanej 


= (jM  - 1.5)  AX 


nM  B 


(4.17) 


XN 


'N 


1 


^u  " (^-  - tan  BB) 

N 


= (jn  “ 1.5)  AX 


Therefore  for  a given  number  of  grid  lines  outboard  of  the  edge,  or 
equivalently,  if  given 


1.5 

175 


(4.18) 


1 


i 

l* 


■ 


the  location  Cu  is 


r _ e t Jmn  ~ 1 , 

“ CR  ' 7 771 

Jmn  nw 


(4.19) 


. 1 ' 


Conversely,  a given  ( £ /?R)  implies  some  jmn>  which  necessarily  should  be 
restricted  to  levels  leading  to  integer  values  of 


jm  - 1.5  + (J„  - l.B)  Jmn 


i.i!.  i 

which  poses  no  difficulty.  Essentially,  eq.  (4.19)  shows  that 


(4.20) 


Jmn  > ^ iMies  ~>Cu>5r  ( - CM) 


(4.21) 


and  the  £u  such  that  jN  is  normal  to  the  upper  (expansion)  surface 
extended  (i.e.  TM)  fixes  x(jM)  = tanBR»  and  from  eq.  (4.17)  is  then 


<g 

u N 


sin3g  cos3g 


In  that  case  also 


x(jM) 


sin  a,  cos 


(4.22) 


(4.23) 


The  compression  side  grid  extends  to  the  intersection  of  the  (i,j)  = (Um) 
meshlines  or  to  (point  P,  Figure  6)  J 


j 


(4.24) 


22. 


’max 


Uu-^)x 

1 + 9’  XM 


M 


The  j line  procedure  may  be  summarized  as  follows: 

a)  Assume  in  order  to  fix  the  mesh  ooint  location  outboard  of 
the  airfoil  edge,  and  from  eq.  (4.22)  the  approximate  (CU)N 


for  a normal  meshline  follows. 


b) 

c) 


Assume  CR/Cj  to  fix  nR  from  eq.  (3.30),  and  also  nM(=  £H/tanBg): 
Assume  jn  and  fix  integer  value  of  j from  eqs.  (4.19),  (4.20). 


The  i line  procedure  is  straightforward  from  eqs.  (4.6)  and  (4.12). 

i , i U ! ■ ' 

tan6 


¥(1) 


(.!n=o  ~ 


nc 


) + «K1)  = (i-1  )Aip  + ^(1) 


(4.25) 


*e(i>  = S=0 


= ( i - 1 ) Aipe  + ipe(l) 


in  which  i^(l)  = -g‘  and  ^e(l)  = 0.  The  choice  of  a maximum  |?n_0l 
based  on  expected  shock  strengths  (see,  e.g.,  eq.  (2.7))  then  fixes 
both  4<(im),  ^g(im)  levels. 

Boundary  Conditions 

Boundary  conditions  are  of  four  kinds:  surface,  symmetry,  overlapping 
computational  grids  and  grid  external  field  edges.  The  latter  were  main- 
tained at  free  stream  property  levels  (along  the  i and  jm  boundaries), 
excepting  for  certain  attached  shock  cases  that  took  advantage  of  the 
known  uniform  regions  of  compression  near  the  edges  to  reduce  consideration 
to  solely  the  lower  region. 


23. 


The  surface  boundary  condition  from  eq.  (3.31)  is 


(— ) - tan 6 

u w 

9' 


(4.26) 


for  the  compression  surface  with  = (tanBp  - tan6)  / s*.  The  same 
form  applies  for  the  expansion  side  with  6 taken  equal  to  zero;  i.e. 


(-) 


w 


tanBg 


(4.27) 


In  practice  such  boundary  conditions  were  enforced  after  each  integration 
step  as  the  calculation  proceeds  and  eq.  (4.26)  indicates  some  freedom  in 
choosing  u,  v,  or  w as  the  basis  for  such  adjustments.  Since  q2  = u2+v2+w2 
(and  omitting  ( ) for  convenience),  eq.  (4.26)  implies  quadratic  equations 

W 

for  v and  w: 


aw2  - 2bw  + c + (q2  - u2)(l  - g12)  = 0 
av2  + 2bg'v  + c = 0 


(4.28) 


Here 


a = 1 + g'2 

b = u tan5 

c = u2(l  + tan26)  -q2 


(4.29) 
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as  is  also  clear  from  the  slopes  of  the  constant  u line  and  its  normal. 

Point  C in  Figure  8,  i.e.  (v,w)  = (0,  utanS),  is  the  reentrant  corner 
condition  as  well  as  the  entire  compression  surface  condition  if  on-design. 

For  off-design  eq.  (4.26)  implies  an  interval  0 < |v|  £ | v T j along  the  ACB 
line  for  some  constant  u.  An  initial  guess  for  the  field  does  not 
necessarily  correspond  to  the  correct  interval,  although  in  the  attached 
shock  case  the  exact  interval  limits  are  predetermined.  However,  in  either 
case  the  algorithm  eqs.  (4.2),  (4.3)  does  not  ensure  consistent  state 
properties  for  eq.  (4.26)  as  the  calculation  proceeds. 

Thus  the  E*  vector  as  given  by  a corrector  algorithm,  eq.  (4.3),  may 
result  in  a (v,w)  pair  that  does  not  lie  on  the  u = constant  locus  (Figure  8), 
and  in  fact  may  lie  within  the  interior  of  the  = 1 circle.  That  the 
possibility  is  of  importance  may  be  realized  by  noting  that  g1  can  be  small, 
so  that  the  points  min  and  C are  in  fact  close  together.  Continuation  of  the 
calculation  requires  distinct  procedures  for  the  k^  < 1 situations.  If 
k^>  1 the  imposed  boundary  condition  after  each  corrector  step  is  given 
by  eq.  (4.33);  if  k <1  the  effective  k may  be  ammended  to  that  satisfying 

i i 1 

eq.  (4.33)  for  v,  namely 

k = 1 + ( V --- -Vin.  )2  (4.35) 

1 wmin 

and  w then  follows  from  eq.  (4.33).  This  procedure  effectively  increases 
w without  constraint  on  the  sign  of  v. 

Alternatively,  with  v as  the  basis  after  each  corrector  step 


- 


- 2vg'tan6  + 'f  (2vg'tan<5)z  - (l+tan26)[(l+gl2)v2-q2] 
1 + tan26 


w = u tan  6 + vg* 


for  the  compression  surface,  and 


V - < 5%  >: 


w = v tan 


i I ! • i 


for  the  expansion  side.  Both  eqs.  (4.33)  and  (4.36),  (4.37)  were  used 
successfully,  the  former  for  the  attached  shock  and  the  latter  for  the 
detached  shock  applications. 

The  symmetry  plane  (q  = 0)  condition  imposes  a reflection  of  the 
j = 1 mesh  line  properties  properly  interpolated  to  take  account  of  the 
diverging  grid  in  the  (n,0  plane.  With  £ and  ip  from  eqs.  (4.7)  and 
(4.25),  and  x(j)  = + Ay/2  for  j = 1.2  the  interpolation  follows  from 
(Figure  9)  the  weighting 


(4.36) 


(4.37) 


= £(1..2)  - ?(i,1) 

€(i  + 1.2) - £(i,2J 


<0  ift^(i)  ^0 


(4.38) 


E.g., 


P(1.1)  = p(i .2)  - [p(i  + 1,2)  - 


P(1.2)]  f§ 


(4.39) 


J 
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Similarly,  p,  u,  w symmetry  and  v asymmetry  may  be  imposed  at  all  (i,l) 
points.  For  the  expansion  side  (£  < 0) 


A?  = e(i,2)  - £(i  + 1,2) 


(4.40) 


p(i ,1 ) = P(i>2)  - [p(i ,2)  - p(i  + 1,2)]  || 


(4.41) 


Outboard  of  the  tip  (n  > 1)  either  symmetry  or  overlapping  grid 
interpolations  were  imposed.  For  attached  shock  conditions  TM  (Figure  6) 

i i i M • • i 

may  be  considered  to  be  a symmetry  plane  and  the  grid  values  for  p,  p,  u, 
and  (v2  + w2)  were  upgraded  after  each  sweep  of  the  entire  mesh  region 
so  as  to  impose  that  symmetry.  For  detached  shock  conditions  the 
overlapping  grid  results  were  used  to  upgrade  the  boundary  conditions 
along  either  TM  or  TP  (for  the  upper  or  lower  region  respectively)  based 
upon  the  results  of  a set  number  of  sweeps  for  the  prior  integration 
region. 

In  the  symmetry  case,  for  each  j _<  j <_  jm-l  the  ip  corresponding 
to  each  x(J)  along  TM  is 


= (SU-CR)  xtangB  - gR 
nR  + x(nRtanp  - Cu) 


(4.42) 


which  fixes  the  maximum  (integer)  i for  a given  j in  region  TPM  on 
use  of  eq.  (4.25).  Since  the  normal  to  TM  has  the  slope  (-  cotBg) 
the  symmetrically  opposite  point  corresponds  to 


(4.43) 


where 


°PP  «u  - «opp 


opp  nR  - nopp 


tan3B  + n.j 

opp  1 + tan3D  tan(2B  - 3.) 

D D I 


f'opp  nopp  tan(26e  " V 


and  3.j  = tan”  (C^/n^ ) - Thus 


P(i.j)  = P(i+l»j)  + [p(i+l»j)-p(i+l»j-l)] 


, , V i.l  ! • • 


X„nn-X(j) 


i +1 ) - \p 

+ [p(i+l ,j)-p(i+2,j)]  2£P 

Aip 


with  similar  relations  for  p(i,j),  u(i,j),  and  vopp,  wQpp.  Then 


(4.44) 


(4.45) 


= vopp  cos2%  + woppsin  2% 


w(i’j)  = vopp  sin26B  “ woppCOS2eB 


such  that  vi 2  + w2  = v 2 + w 2 . The  levels  from  eqs.  (4.45)  ,(4.46) 

opp  opp  ^ 

are  then  associated  with  all  i equal  or  less  than  that  from  (4.42)  to 
complete  the  grid  specification  throughout  TPM. 


(4.46) 
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The  procedure  for  the  overlapping  grid  case  is  similar  to  the  above. 
Partially  converged  results  for  the  lower  region,  e.g.,  provide  conditions 
along  TM  as  the  tip  condition  for  the  upper  region  integration.  Eq.  (4.42) 
defines  both  for  a given  jn  < j < j -1 , and  with  eq.  (4.25)  the 
adjacent  i values  (i+,  i~)  bracketing  the  point  on  the  TM  locus.  Then 


pTM  •=  Pj-  + (p,+  - PrM 


(4.47) 


and  similarly  for  the  other  state  properties.  On  completion  of  partial 
convergence  for  the  upper  region,  the  same  procedure  establishes  values 
along  TP  as  the  subsequent  lower  region  condition.  In  that  case  the  \p£ 
corresponding  to  each  x(j)  point  on  TP  is 

tan6(l  + x[tan3B  - £,,]) 

= (4.48) 

6 s*(l  + g'x) 


Decoding 

The  flow  variables  must  be  decoded  from  E*  or  E*  results  in  order 

e 

to  reconstruct  F*,  G*,  H*  after  each  predictior  and  corrector  step.  From 
eqs.  (3.2)  and  (3.4) 

E + ^ E2  - 4k(l-k) (E2  - E2  - E2) 
u = _2 2 1 3 *» 

2(l-k)E1 


v = 


E 

3l 

E 

l 


(4.49) 
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(4.49) 

E (continued) 

W - £i 
1 

Ei 

P = T 

P.:  = ^ 0-u2  - (f-  )2  - (^)2  ) 

i i 


The  sign  choice  for  the  radical  follows  from  consideration  of  small  v 
and  w,  or  that  u remains  supersonic.  Once  again,  however,  the  finite 
difference  procedure  does  not  ensure  consistent  E.  values  after  each 
Ac  step.  If  R is  the  contents  of  the  radical  in  eq.  (4.49)  there  are 
two  implied  physical  constraints: 


R > 0 

E + AC 

< 1 

2(l-k)Ei  “ 

i 

The  first  ensures  real  u and  the  second  that  u _<  q _<  1.0. 
I.e. , 

= (4k(l-k)[l  - (q2-u2)])1/2 

min 

= 1 - k(qz  - u2) 
max 


is  the  permissible  E^/E^ range  for  the  calculation  to  continue. 

< 


(4.50) 


(4.51) 


. . .. 
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These  bounds  are  of  special  importance  for  undershoot/overshoot  behavior 

near  shock  locations,  for  near  sonic  speed,  and  at  high  Mach  numbers. 

In  the  present  computations  both  the  R < 0 and  p < 0 possibilities  were 

monitored  and  corrected  for  by  modification  of  E to  constrain  the  system 

2 

, * i 

to  physically  meaningful  results  during  the  asymptotic  approach.  Such 

modifications  proved  necessary  when  starting  with  arbitrary  initial 

conditions  but  the  need  does  decay  (i.e.  in  terms  of  number  of  necessary 

modifications)  as  the  integration  proceeds. 

Figure  10  illustrates  the  limits  of  eqs.  (4.51)  and  also  includes 

the  E /E  variation  with  M for  the  special  case  of ■ v = w = 0;  i.e.  in 
2 1 

the  latter  case 


[IiIm2  (i+1i1m2)]1/2 


(v  = w = 0) 


(4.52) 


1 / 2 

and  varies  between  ( y2- 1 ) /y  = 0.6999  and  unity  for  1 < M <_  In 
general  the  minimum  condition  corresponds  to 


u 2 > l^-1-:-q2-)  (4.53) 

~ 1 -2k 

and  the  v = w = 0 special  case  implies  u = / k/(l-k)  or  sonic  speed, 
consistent  with  eq.  (4.52). 


Solution  Procedures 

The  integration  of  the  field  proceeds  from  assumed  initial  conditions 
at  the  arbitrary  axial  plane  c = 1.  Free  stream  conditions  suffice  but  an 


appreciable  reduction  in  computation  time  results  with  either  imposed 
shock  layer  flow  variables  which  are  based  on  psuedo  on-design  conditions 
or  the  continuation  from  a previously  determined  solution  as  the  conditions 
closely  approximating  those  next  of  interest.  An  attached  shock  evaluation 
would  therefore  typically  follow  by  taking  a aM^  increment  from  an  already 
converged  solution  for  some  Mro,  starting  originally  with  the  on-design 
case.  At  each  c the  predictor  and  corrector  computations  from  eqs.  (4,2), 
(4.3)  imply  decoded  values  applicable  at  c + Ac.  boundary  conditions  are 
then  enforced,  and  the  steps  are  repeated  until  uniform  conical  conditions 
are  achieved  at  some  downstream  t,. 

. , ill!.  i 

For  known  attached  shock  conditions  (or.Cr)  was  placed  at  the  side 
edge  and  £u  moved  to  00  such  that  the  j mesh  lines  formed  a set  parallel 
to  the  n = 0 centerline. In  such  cases  only  the  lower  (compression)  region 
need  be  considered  and  the  outermost  grid  boundary  (jm)  was  displaced 
inboard  so  as  to  lie  within  the  uniform  flow  region  between  the  tip  and 
the  post  shock  Mach  cone  from  the  apex  (Figure  11).  This  arrangement 
provides  an  i grid  line  which  parallels  the  shock  in  the  uniform  region 
and  places  the  singular  ( )^  origin  external  to  the  computational  domain. 
However,  free  stream  conditions  are  no  longer  appropriate  along  jm  and  were 
replaced  in  such  cases  with  the  simple  floating  uniformity  condition 
p(i.jm)  = p(i.Jm-l).  etc.  The  tip  region  properties  therefore  were  self 
adjusting  to  a final  sweptback  edge  compression  level. 

For  detached  shock  conditions,  as  well  as  when  check  computations 
of  the  attached  shock  case  were  considered  the  alternating  grid  procedure 
was  employed  with  partial  convergence  of  the  respective  upper  and  lower 
regions  carried  out  successively.  Initial  conditions  consisted  of  finite 
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shock  layer  approximations  extended  over  the  entire  lower  (compression) 
region  and  the  tip  portion  of  the  upper  region.  This  is  of  some  importance 
for  small  6 since  the  overlap  region  TMP  in  that  case  may  itself  be  small. 


I 
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5.  NUMERICAL  RESULTS 

Computations  were  carried  out  for  two  basic  caret  sections  with 
= 35°  and  6 = 5.746°  and  22.181°  (Table  3).  The  geometries  correspond 
to  on-design  flows  at  Mg  = 2 and  4 respectively  and  in  each  case  both 
s*  = 1 and  2 were  considered.  For  the  Mg  = 2 geometry  both  the  attached 
and  detached  (alternating)  algorithms  were  used  for  some  of  the  attached 
side  edge  cases  as  a consistency  check  on  the  procedures. 

A measure  of  the  shock  capturing  behavior  is  shown  in  Figures  12  a, b 
for  several  < Mg  at  each  of  the  Mg  levels.  The  pressure  distribution 
through  the  shock  layer  adjacent  to  the  symmetry  plane  (n  = 0.01)  for 

, , ill*.  1 

attached  side  edge  cases  exhibits  fairly  uniform  conditions  both  within 
and  external  to  the  shock  layer.  Overshoot  and  undershoot  was  always 
present  upstream  and  downstream  of  the  shock  location  with  a numerical 
adjustment  to  the  pressure  jump  occurring  over  3 to  4 mesh  widths.  The 
major  shock  jump  was  over  two  mesh  widths  as  indicated  by  the  dashed  lines. 

Lower  (compression)  surface  pressure  coefficient  distributions  appear 
in  Figures  13  a,  b,  c.  Figures  13  a and  c contain  only  attached  side 
edge  cases;  Figure  13  b includes  those  for  side  edge  detachment.  The 
attached  algorithm  imposed  the  uniformity  condition  at  n > n$  (Figure  11) 
and  led  to  converged  solutions  in  good  agreement  with  the  theoretical 
pressure  levels  (indicated  by  horizontal  lines  in  Figures  13)  for  the 
uniform  edge  regions,  and  with  break  points  in  the  distribution  quite  con- 
sistent with  the  expected  ns  locations.  The  detached  shock  cases 
(M^  = 1.3,  1.4  in  Figure  13  b)  achieved  near  edge  pressures  comparable 
with  post  shock  levels  for  a symmetric  shock  lying  at  1 < n < 1.06. 
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E.g.,  for  the  M =1.3  case  0.38  < C„  < 0.41  is  implied  behind  such 

J 00  P r 

shocks  when  lying  outboard  of  the  side  edge. 

Upper  (expansion)  surface  pressures  for  detached  cases  indicated 
free  stream  levels  over  the  entire  surface  span.  The  expansion  about 
the  side  edge  rapidly  weakens  the  nearby  shock  and  resulted  in  a pressure 
field  for  n < 1 within  approximately  one  percent  of  the  free  stream 
pressure.  Nevertheless,  the  upper  region  is  of  some  importance  as  a 
boundary  in  the  detached  region  which  allows  for  the  rapid  changes. 

Shock  jump  locii  in  the  crossplane  are  shown  in  Figures  14  a,  b,  c, 
the  geomet'ry  being  inverted.  All  points  were  obtained  as  midpoint 
estimates  of  the  greatest  pressure  jump  location  along  each  grid  line 
(e.g.,  as  in  Figure  12).  Solid  lines  in  Figures  14  a,  c are  the  theo- 
retical attached  shock  traces  for  n > ns  (Figure  11).  The  interpolated 
points  are  in  good  agreement  with  such  directions  and  maintain  the  same 
direction  for  an  additional  interval  n < n$  in  which  the  edge  conditions 
persist  away  from  the  surface  and  exterior  to  the  n$  Mach  cone. 

The  dashed  line  in  Figures  14  a,  c indicates  the  limiting  shock 
directions  for  post  shock  supersonic  flow  normal  to  the  edge.  For  these 
geometries  (Table  3)  the  condition  corresponds  to  = 1.428  and  2.326 
for  the  (6,  6g)  pair  examples.  The  = 1.4  case  in  Figure  14  b represents 
a condition  very  close  to  the  classical  detachment  level  (M^  = 1.41); 

M^  = 1.3  is  a fully  detached  case  from  the  side  edge  and  is  approaching 
close  to  the  apex  detachment  limit  of  M =1.26. 

In  both  Figures  14  a and  c the  point  scatter  increases  appreciably 
as  the  shock  traces  depart  from  the  computational  grid  directions, 

) 
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i.e.  for  n -*■  0 regions  and  for  increasing  |MD  - M^J , with  an  indicated 
waviness  on  the  order  of  the  grid  scale.  For  the  detached  M =1.3 
case  the  grid  scale  was  approximately  halved  to  confirm  the  effect  for 
a widely  bulging  shock. 

The  Mro  = 1.4  case  closely  adhered  to  the  limiting  shock  direction 
for  a supersonic  side  edge  with  virtually  no  discernible  pressure  jump 
in  the  field  on  the  expansion  side.  For  M^  = 1.3  the  detachment  is 
quite  clear  and  an  appreciably  weakened  shock  is  visible  on  the  expansion 
side.  Two  Mach  cone  traces  are  shown  for  the  upper  (expansion)  surface 

in  Figure  14  b.  One  is  the  apex  Mach  cone  for  M =1.3  which  extends  to 

• . « » * ! * « 

about  n = 0.6.  The  other  is  the  free  stream  Mach  cone  disturbance 
envelope  limit  from  the  side  edge  which  is  tangent  to  the  Mach  cone 
centered  on  the  apex.  The  latter  limit  makes  clear  the  rapid  weakening 
of  the  shock  in  the  close  vicinity  of  the  edge  and  the  reasonably 
parallel  shock  trace  extending  into  the  upper  field.  This  appears  to 
be  consistent  with  the  virtually  undisturbed  surface  pressures  found 
for  that  surface  as  mentioned  above.  The  detachment  scale  itself 
remains  small  (n  ~ 1.05)  despite  the  close  approach  to  complete  detach- 
ment at  M =1.26. 

OO 

Figures  15,  16,  and  17  summarize  the  overall  force  coefficients  for 
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with  (Pbase/P<J  taken  to  be  zero.  The  base  pressure  level  can  make  a 
substantial  contribution  to  a drag  reduction,  of  course,  depending  on  the 
Mach  number  level.  For  example,  on-design  the  extremes  of  ( P^/p^)  for 
Md  = (2,4)  imply  the  range  CD  = (0.0492,  0.2157)  to  (0.0133,  .1793). 

For  M = 2,  C, /Cn  increases  from  2.68  to  9.94,  and  for  M = 4 from 
2.04  to  2.45  for  (p^/p^)  in  the  range  0.  to  1 . 

Comparison  is  made  between  results  from  the  attached  and  detached 
algorithms  in  Figures  15  and  16  for  s*  = 2 and  M^  = 1.7  and  1.5.  The 
disagreement  for  the  latter  may  be  attributed  to  edge  effects  on  the 
upper  surface  under  near  detachment  conditions,  as  evident  in  the  Cn 

, , ; 1 1 ! • i u 

unaffected  agreement  for  the  same  cases.  For  relatively  smaller  n$ 

(M^  = 1.7)  the  alternating  procedure  proved  to  be  very  consistent. 

The  coefficient  variation  in  the  detachment  range  of  Mro  shows  no 

abrupt  lift  losses  due  to  edge  effects.  One  computation  was  purposely 

carried  out  beyond  the  apex  detachment  limit  as  a consistency  check; 

for  Hm  = 1.2  a sharp  decrease  in  CL  was  indicated  (approximately  20%) 

on  starting  from  the  M^  = 1.3  solution  and  completing  several  hundred 

♦ 

Ac  steps.  The  lift  to  drag  ratio  decreases  smoothly  from  2.68  to  2.48 
for  attached  shocks  at  2.0  > Mro  >_  1.50  and  averages  2.47  ± 1%  in  the 
detached  M interval. 

oo 

The  detachment  evaluations  depend  critically  on  an  appropriate 
initial  condition  in  the  overlapping  mesh  regions  PTM,  Figure  6.  For 
thin  geometries,  such  as  the  6 = 5.75°,  MD  = 2 example  cited  here,  it 
is  essential  that  initial  conditions  along  the  boundaries  PT  and  TM  be 
representative  of  a shock  layer,  and  that  the  overlap  be  sufficient  to 
I 

i 
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allow  description  of  the  gradients  in  that  region.  For  P M,  i.e.  an 
extremely  thin  geometry,  the  alternating  algorithm  would  effectively 
require  a match  of  derivatives  along  the  common  interface.  In  this 
sense  the  calculations  for  6 ~ 5°  are  a good  test  of  the  capabilities  of 
overlapping  mesh  regions  chosen  so  as  to  conform  to  specific  surfaces 
and  shock  layer  segments.  It  is  to  be  expected  that  "thick"  caret 
solutions,  e.g.,  would  converge  somewhat  faster  and  be  less  sensitive 
to  initial  field  choices.  Virtually  all  of  the  present  results  were 
obtained,  however,  by  stepping  to  a new  level  and  proceeding  from 
the  previous  field.  , 
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6.  CONCLUSION 

The  study  has  applied  shock  capturing  finite  difference  procedures 
to  waverider  configurations  of  the  caret  type.  It  is  fairly  evident  that 
a grid  definition  for  crossplane  sections  with  cornersand  concave  regions 
poses  special  constraints.  For  off-design  Mach  numbers  such  that  shocks 
remain  attached  at  side  edges  it  has  been  shown  that  computational  grids 
may  be  chosen  to  conform  to  crossplane  surface  and  symmetry  planes  while 
simultaneously  providing  approximate  but  adequate  matching  with  antici- 
pated shock  traces.  For  the  Mach  number  interval  corresponding  to  edge 
detached  to  apex  detached  flow  fields  the  interacting  compression  and 

• i i h !•  ' 

expansion  portions  of  the  field  may  be  evaluated  by  similar  specialized 
grid  regions  of  relatively  small  overlap. 

The  suggested  alternating  algorithm  requires  sequencing  the 
separate  integrations  for  each  such  region,  with  interpolative  adjustment 
of  the  boundary  conditions  in  the  edge  region  as  the  asymptotic  integration 
proceeds.  The  evaluations  demonstrate  the  utility  of  such  an  alternating 
algorithm  for  relatively  small  regions  of  overlap  that  are  implied  by  thin 
configurations. 

While  no  attempt  was  made  to  refine  grid  scales  so  as  to  define  the 
flow  details  in  the  detached  near  edge  region,  the  overall  force  coefficient 
results  make  clear  that  sharp  discontinuities  in  behavior  are  not  present 
in  the  transition  from  attached  to  detached  side  edge  conditions. 
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APPENDIX  1 
COMPUTER  CODE 


A code  glossary,  flow  chart  and  listing  are  given  below  for  the 
detached  shock  field  algorithm.  Computations  were  carried  out  on  the 
M.I.T.  Multiplexed  information  Computing  Service  (MULTICS)  system 
which  makes  use  of  a Honeywell  6180  computer  in  a time  sharing  mode. 
Computer  time  per  grid  point  was  2.34  x 10_3seconds  and  solution  time 
varied  from  10  minutes  for  nearly  on-design  to  two  hours  for  a fully 
detached  case. 

,,  4 u !•  i 

Specific  grid  constants  are  summarized  in  Table  4 for  the  cases 
considered.  Stable  solutions  were  completed  with  bet  = 0.05. 
Numerical  instabilities  were  evident  for  bet  > 0.10. 


. . 
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Code  Glossary 

a 

(1  + tan26) 

ajmn 

(VT.Sj/Un-l.S) 

ak 

(Y-D/2Y 

akl 

nM 

ak4 

initial  condition 

al 

Ac/ Ay 

am 

M 

amd 

oo 

on  design 

b 

vh  tan6 

be 

3D 

bet 

Ac/  At/; 

chi 

X 

d 

(l+h2)v2-q2 

dz 

dC 

dec 

AX 

dep 

Ai/; 

del 

6 

del  chib 

x(jm)  - xO) 

del  psib 

^(im)  - ¥(1) 

dxir 

6S/AS 

e 

e • ; 

et 

n 1 

etb 

0 along  0B  locus 

etl 

\ 

etr 

nR 

f 

F 

g 

G 

ga 

Y 

h 

H 

hh 

(tan&D  - tan6)/s* 

im,  iml 

i 

iii 

m 

1-1 

V j"1 

m 

jm 
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) 


jn 

jn 

jnl 

jn  + 1 

jjj 

Jm  " 1 

ni 

cycle  count 

P 

P 

Pi 

7r 

PP 

(P/PoL 

psi 

Pi 

P/P0 

P2 

P/P»" 

q 

q 

qq 

q2 

r 

p 

rr 

(p/p0L 

rl 

p/p0 

s 

4-  °° 

s* 

tbb 

tangg 

tbe 

tanBD 

tbbi 

tan(2Bg  - ) 

u 

u 

uu 

V 

V 

w 

w 

xi 

5 

xib 

£ along  Bg  loc 

xir 

?R*  ?M 

xiu 

Su 

xid 

Sim  for  ^ = 0 

z 

C 

dimension  p(42,42),r(42,42),u(42,42>  , v ( 42 , 42 ) , w ( 42 , 42 ) 

dimension  et (42,42),  xi  (42,42) , ch  i (42) , psi  (42) 

dimension  pi  (42, 42), rl (42,42)  , u I (m2  ,42 ) , v I (42, 42)»wl (42,42) 

dimension  e 1 1 ( 42  » 42) »xil(42,42)»chil ( 42 ) * ps i I (42) 

dimension  e(42,42,4),f(42,42»4),g(42,42,4)»h(42,42,4) 

dimension  el(42,42,4) 

common  /cl/  p,r,u,v,w 

common  /c2/  ga , am , be , p i , amd 

common  /c3/  hh , dec , dep , tde, occ 1 , dep I 

common  /c4/  i m , J m , im I , ) m I 

common  /c5/  ak9,ak4 

common  /c6/  ] ) n, a, xi u, tbb, etr  i 

common  /c7/  i i i , J ) ) , i i i I , ) ) 1 I , 1 nl , | t , ] 1 1 
common  /c8/  xi  , et,xil,etl 
common  /c9/  )n 

common  /clO/  xi r , chi , psi , chi  I , ps i I 

common^/cll/  ak 

common  /cl3/  s 

common  /cl  4/  el,e,f,g,h 

common  /cl5/  dz,al,bet,z 

common  /ci6/  del chip,  delpsib,del chi  I 

common  /cl  7/  neqr,ne3p 

common  /cl  8/  p 1 , r I « u I , v I , w I ■;  fu''  1 

common  /cl9/  nit 
common  /c2  0/  j a , ] b 

namelist  /nan/  ga,am,StDe,del ,bet, |nt etal , 
akl,ak2»ni,amd,ak3,nit,ak4 

rewind  22 

constants,  coordinates 


read 

( 22,nam) 

ak 

= 

( ga-1 . ) / ( 2,  *g  a ) 

Pi 

= 

4 .*  at  an  ( 1 . ) 

tde 

= 

tan  (del/180.*pi) 

a 

— 

1 tde**2 

tbe 

= 

t sn (d  e/18  0 * *p  i ) 

tbb 

= 

tte/s 

hh 

= 

(tbe-tde)/s 

xir 

= 

akl * tbb 

etr 

= 

(xi r - ( tde/s ) ) /hh 

x i u 

2.*etal/(sin(2.*atan(tbb)  )) 

a J mn 

( (xiu/etal ) -tob) /( (xid/akl)-tbb) 

ak9 

1 .5  + ( J n-1.5  > *ajmn 

1 m 

i (i  x ( ak9 ) 

i m 

= 

J ir 

] m 1 

- 

) t 

i m 1 

-= 

i IT 

a ) mn 

— 

( Jm-1.5) / ( J n-1.5) 

xiu 

xir*(a]mn-l.) /(a]mn-(akl/etal  ) ) 

x i c 1 

ak2  *x  ir 

dec 

l./((xiu/akl-tbb)MJm-1.5)) 

dep 

(xicl-(tde/s)  )/etr/(im-l*  ) 

dec  1 

= 

dec 

dep  1 

ak3/( iml-1.  ) 

] 1 J 1 

= 

j rrl  - 1 

iii  1 

ini  - 1 

ini 

-r 

] r ♦ 1 

4b 


] 1 n = 1 r-  1 
dz  = bet/iii 

al  = d i*  J J ] 

print  10C»  am, d e I ,be, ( lft 0 . /p i ) * at an ( tbb ) * s • 

& e t a I » akl  , e t r , x iu  , xi  c I , 

l ak2 , dec, Jcp  , ak4 , 

& ak3, dccl , dcp I « 

1 Jn. J m, im, 

K ni t , ) m I , i m I , 

1 al »bet,amd 

c 

doll  i = 1 , i m 

psi(i)  = (i-1.  )*dcp  - hh 

do  10  | = li|n  ; ' 1 

chilj)  = (J-1.5)*dcc 

ak9  = 1.  - chi ( | ) *psi ( i ) 

et(i,J>  = chi ( J ) * (x i u-xi r-e t r*psi ( i ) ) /ak9 

10  xi  ( i , 1 ) = (xir  + (etr-xiu*cni()))*psi(i))/ak9 

11  continue 

do  204  J = 10,IJI 
ak9  = etCl, J ) - 1. 

J tl  =1 

if  (ak9.gt.O.)  go  to  205 

204  continue  , 

205  J t = Jtl  - 1 : ; 

print  146,  J t 

print  120 

delchib  = chi(Jm)  - chi(l) 

delpsib  = psi(in)  - psi(l) 

do  61  i I = i, i m I 
psil(il)  = (il-l»)*dcpl 

do  62  ) I = 1, J ml 

chi  I ( J I ) = ( J 1-1.5) * dec  I 

ak9  = 1.  + chil(JI)*tbb 

etl (il  ,J  l)=  cnil  (]  I ) * ( xiu-psi 1 (i  I) ) /ak9 
62  xil(i!,JI)=  (chi  I ( J I ) *xiu*tbb  + psi I ( i I ) ) /ak9 

61  continue 

delchil  = chil(jml)  - child) 

call  prcoord 

c initial  conditions!  -0+  disc,des,-n 

print  102 
read  108,1 
call  ic  ( I) 

c integral  ion 

z = 0 . 
kc  =0 

208  print  147 

read  108,1 

23  do  21  k2  = l,ni 

kc  = kc  ♦ 1 

J a = 1 1 

lb  = I t - 1 

do  207  kl  = 1 , ni  t 

call  int  eg  rated) 

207  continue 

la  = 2 

Jb  = 1 ’ 

call  integrated) 

21  continue  ' 

print  120  . 


J 


24 

26 

25 

49 

51 

50 
c 

100 


102 

105 

106 

107 

108 
109 
120 
135 
124 

146 

147 


pri n i j.  <•  ■> $ ftC i i 

print  124,  negr, negp 
negr  = C 
negp  = C 
z = z + dz*ni 

not  converged:  -0+  mane yc « st oo , se I mod « -0+  orf  intstop 

rewind  9 

write  (9)  p , r, u, v , w, o I , r I , u I , v I , wl 

print  106iZ 

call  pr 

print  107 

read  108,1  I 

if  (in  23,25,26 

cal  I modd  ) , i 

go  to  208 

print  109 

read  108,1  I 

if  (II)  49,50,51 

print  135 

call  pr  f in 

continue 

stop 


& 

S> 

k 

k 

k 

k 


f 6 .2 

•• 

t 

de 1 tu=", 

• 

CO 

4," 

b* 

- 

", f 5.1 ,/, 

"et 

cl  - 

— •» 

— t 

t 8 . 

4,"  xiu 

•• 

,f  8 

, f 6. 

2 , " 

delchi-" 

,f  8 

• 4y 

, f 6 . 

2," 

de 1 c h 1 = " 

,f  8 

• 4 t 

, i4. 

4 X , 

" J m =", 

i 4 , 

6 x * 

, i 4,  4x  , 

“ J ml  =”, 

i4, 

fe  X t 

t f 6 • 

2," 

oet 

, f 8 

• *4  t 

format  ("m  =' 

f 8.4,' 


" ] n 

"nit  =* 

"a  I =* 
f ormat  ("fix  i .c .") 
format  (*’i  ters="  » i 7,  i 7) 
format  ("In  x =",f8.2) 
format  ("proceed  option") 
format  ( i2  ) 

format  ("terminate  option")  . 
format  (2x1 

format  (/,/,/)  ; 

format  ("negr  -" , i 6 , 2 x, "negp  ="»i6) 
f ormat  ("J  t =",  17) 
format  ("fix  region") 


be t ad=" , f 8 • 4, " betab- 
(6k  2 * • etam  =“,  f 8.4, 

> xicl  ="  ,f  8.4,/, 
del psi=", f8 .4,"  aK4 
delpsl=",f8.4,/, 
i m ="  , i4 ,/ , 
iml  =",i4,/« 
amd  =" ,fb.4,/) 


c 

c 

c 

c 


end 


subrout i re  ic ( I ) 

initial  conditions 

dimension  p(42,42),r(42,42) ,u(42,42) , v ( 42 , 42) , w (42, 42 ) 

*—  "* (42, 42),rl (42,42) ,u! (42,42) ,vl (42,42 >,w I (42,42) 

p,r,u,v,w 
ga,am,be , pi , amd 
hh,dcc,dcp, tde,dccl ,dcpl 
im,  ) m,  im  I , J ml 

aK9,aK4  i 

iii, J | |, i li  I, ) n I ♦ jnl, ) t, j tl 

] n I 

pl,rl,ul,vl,wl 


dimension  pi 
common  /cl/ 
common 
common 
common 
common 
common 
common 
common 


/c2  / 
/c3  / 
/c4  / 
/c5/ 
/c  7/ 
/c9  / 
/cl  8/ 


, wr  i t e 


t 


j , C 1 . i.  1 ' V <?  1 s 

rr  = (1 4 ♦ ( y a-1.  ) * am*  am/2  « ) * * (1  • / ( 1 . -ga  ) I 

pp  = rr**ga 

uu  = sqrt  ( ( ga- 1 . ) ♦am* a m/2 . / { 1 . 4 ( ga-1 . ) * am *am/ 2 . ) ) 

if  (I)  12,13,14 

c read  disc  field 

12  rewind  9 

read  (9)  p,r»u,v,w,pl ,rl,ul , v I * w I 

call  pr 

return 

c “on  design"  i.c. 

13  i 9 = ifix  (l.+hh/dcp) 

ilO  = ]n  ♦ 3 

do  15  i = him 

do  15  | = 1 , J m 

p ( i , J ) = pp 

r ( 1 , j ) - rr 

u ( i , J > - uu 

v ( i , J ) - 0 • 

15  w(i,  | ) = 0 . 

do  75  il  = 1 , i m t 
do  75  J I = 1 , 1 ml 
p 1 ( i I , J I ) = DU 
r I ( i I , | I ) = rr 

u I ( i I , J I ) = uu  , , t , I 

v I ( i I , J I)  = 0. 

75  wt  ( i I , J I ) = 0. 

ak9  = ( am* s i n ( ( be/ 180 . * p i ) + ( asi n ( 1 . / am) -a s i n ( 1 . /amd ) ) *ak4 ) ) **2 

pi  = pp* ( 1 . *2 • * ga* ( ak9- 1 • ) / (ya +1 • ) ) 

rl  - rr * ( ga  + 1 • ) * ak9/ ( ak 9 * ( ga- 1 . ) +2 • ) 

q = sqr t ( 1 .-pi /r 1 ) 

ul  - q*sqr t ( 1 ./ ( 1 . + t oe* *2) ) 

vl  = 0 . 

wl  = ul*tde 

do  17  i = 1, i9 

do  17  1 = 1 , i 1 0 

pti » 1 ) = Pi 

r < i , J ) = rl  • 

U ( i , J ) = U 1 

v ( i , ] ) = vl 

17  w ( i , J ) - wl 

do  206  i = 1,8 
do  206  1 = J tl, ilO 
pi ( i , J ) = ( pp+pl) /2 . 

r 1 ( i , J ) = ( rr  + rl ) /2  . 

ull(i,|)  = sqrtt  (1. -pi  (i,J)/rl  (i,J  ))/(!.  + tde**2)l 
v*  ( i ,1  ) = 0 . 

206  wj (i,J  ) = ul (i,J ) *t  de 

cal  I pr 
return 

c new  m inf  ini  ty 

14  print  103 
read  104,  am 
rewind  9 

read  (9)  p ,r ,u, v , w,p I , r I , u I , v J , w I 

rr  = ( 1 . 4 ( ga-l. ) *am*am/2 • ) * * ( 1 . / ( 1 . -ga> ) 

pp  = rr**ga  ! 

uu  = sqr  t ( ( qa-1 . ) * am*am/2 • / C 1 . +( ga- 1 • ) * a m* am/2 # ) ) 

do  19  I : 1,1m  , 

p ( 1 m , 1 ) = pp 


r ( i m , J ) = 

rr 

19 

u C i m , J ) = 

uu 

do  76  I 1 = 

1,  J ml 

p 1 ( i m 1 , ] 1 ) 

= OP 

r 1 ( iml  , J 1) 

= rr 

76 

ul (iml  , J 1) 

= uu 

ao  20  i = 

l,im 

p ( i , J m ) = 

PP 

r ( i , J m ) = 

rr 

20 

u( i , | m ) = 

uu 

do  77  il  = 

1,1ml 

pi  ( i 1 , 1 m 1) 

= PP 

r(  ( i 1 , J m 1) 

= rr 

77 

ul  ( i 1 , J ml ) 
cal  1 pr 
return 

= uu 

c 

103 

format  ("m 

= ") 

104 

format  ( f5 

.3) 

c 

end 

c 

c 


88 


c 

c 


79 

c 


ii 


subroutine  integrate(l) 


cal  I 
cal  I 
call 
call 
call 
cal  I 
cal  I 
call 
return 
end 


bell) 
ee  ( I ) 
f gh(l  ) 
prd 

decod  el ( I ) 
f gb  (I  ) 
ccr 

decod  e ( I ) 


subrout  ine  be ( I ) 

boundary  conditions 

dimension  p(42*42)»r(42,42),u(42»42)  , v ( 42  » 42 ) , w ( 42 * 42  ) 
dimension  pi ( 42 , 42)  » r I (42,42) ,ul (42,42) ,vl (42,42),wl (42 
dimension  etl (42,h2I , xi I (42,42) ,chil ( 42 ) , ps i I (42) 
dimension  et(42,42),xi(42,42),chi(42) , psi ( 42) 
common  /cl/  p,r,u,v/,w 

/c3/  hh , dec, dep* tde, dccl , dco I 
/c4 / im, J m, im I , J m I 

)J  n, a , xi  u * t bb , e tr 
ii  i,  I ] 1,  i ii  I , ) J J I , J nl, J t , ] tl 
xi,et,xil,etl 
Jn 

xir,chi,psi,chil,psil 


42) 


common 
common 
common 
common 
common 
common 
common 
common 
if 


/c6  / 
/cl  / 
/c8  / 
/c9  / 
/cl  0/ 


/cl  8/  pl,rl,ul,vl,hl 
(l.gt.O)  go  to  76 
fl  at  surface 


i = 1 
do  7 9 
dq  = 1 . 
ul  ( i , J ) 
wl  ( i,  | ) 


) = 


2, It 

- pi  (i  , ) )/rl  ( i,  J ) i 

= sqrt(qq-(vl (i,| )/cos(atan(tbb) ) ) * *2 ) 
= tbb*vl  (i, J ) i 

centerline  , 

do  80  i = 1,  iii  I 


/•o 

oxir-  - (*ii(i,i)-xil(i,i))/(xil(i,2)-xil(i-ti,2)} 

pl(i,l>  = pi <i,2)+(pl (i+l,2)-ul (i,2) >*dxir 
rl(i,l)  = r I ( i i 2 ) «■ (r I ( i +1 ,2) -r I ( i , 2 ) ) * dxir 
ul(i,l>  = ul  (i,2)+(ul  (i+l,2)-ul  (1,2)  )*dxir 
vl(itl)  =-v/l(i,2)-(vl(i+l,2)-vl<i,2))*dxir 
80  wl(i,l)  = m!  U,2U(mI  (i+i,2) -rtl  (i,2)  ) *dxir 

return 
c 
c 

78  I 

do  30  ] 
qq 

b 
d 

u(ii]  ) 

30  w ( i i 1 ) 
c 

do  31  i 
dxir 
P ( i , 1 > 
r(i ,1) 

Uli,i> 
v 1 1 , 1 ) 

31  w(i,l) 
return 
end 

c 
c 

subroutine  decodel(l) 

dimension  p(42*42)*r(42,42),u(42,42),v(42,42)  , w ( 42  * 42  ) 
dimension  pi  ( 4 2 , 42 ) , r ( (42,42)  »ul  (42,  42 ) ■»  u I (42,42)  ,wl  (42,42) 
dimension  e(42,42,4)  , f (42,42,4)  ,g(>»2,42,4)  ,h(42,42,  4) 
dimension  el(42,42,4) 
common  /cl/  p,r,u,w,w 

common  /c 7/  iii*J  ] J,  iiil, J|)  l,)nl,)t,Jtl 

common  /ell/  ak 

common  /cl4/  el,e,f,q,h 

common  /cl  7/  neqr,ne jp 

common  /cl8/  p I , r I , u I , v I , w I 

common  /c20/  ]a,]b 

akll  - 1 . - ak 

do  36  i = 1 , ii i 

do  36  I = Ja,  ] 1 ) 

if  ( I .gt  .0  > go  to  92 

vl  ( i , J ) = el (i , J ,3)/el(i, J ,1) 

w«(i,J)  = el ( i , J , 4) / el ( i , J , 1 ) 

38  ak9  = el(i,J ,2)/el(i,J ,1) 

aklO  = 4 . *akii *ak* ( 1 . -v I (i , J ) ** 2-w I ( i , J ) * *2  ) 
t = sK9**2  - aklO 

if  (t.ge.O,)  go  to  37 
negr  = negr  ♦ 1 

e 1 ( i , ) , 2 ) = 1. 01*el( i , J ,l)*sqrt ( aklO) 
go  to  38 

37  ul(i,J)  = (ak9+sqrt(t))*0.5/akll 

r 1 ( i , ] ) = el (i , | ,1) /ul  (i,J  ) 

pH  ( i , | ) = rl(l,|)*(l.-ul(i»JI**Z-vl  (i,  J )**2-*il  (i,  1 » **  2 ) 

t = p I ( i , ) ) 

if  ( t. gt.O  * ) go  to36 
negp  = nego  ♦ 1 

go  to  36 


compression  surface 
= 1 

= 2 * J t 

= 1,  - p(i,|)/r(i,j) 

= hh*  t d e* v ( i , ] ) 

= (l.+hh**2) * (v(i,J )**2)  - qg 
= C-b+sqrt (o **2-a*d) ) /a 
= tde*u ( i , 1 ) ♦ hh*v ( i , J ) 
center! ine 
=1 , i i i 

= (xi(i,2)-xi(i,l))/(xi(i+l,2)-xi(i,2)) 
= p(i,2)  - ( o (i +1 ,2) -p ( i ,2 ) ) ♦dxir 
= r(i,2)  - (r ( i +1 ,2 ) -r ( i , 2 ) ) ’ox  i r 
= u ( i , 2 ) - (u(i+l,2) -u(i»2))*dxir 
-- v ( i , 2 ) ♦ Cv (i+1,2) -v(i ,2) ) *dxir 
= w ( i , 2 ) - (w(i+l,2)~w(i,2))*dxir, 


J 


r/> 

v(i*])  = el(i»)*3)/el(i*]«l) 
wfi*J)  = e 1 ( i t J » 4 > /el ( i * J « 1 ) 
ak9  - e 1 ( i i ] ,2)/e±(i,l  il) 

aklO  = 4.*akll*akMl.-v(i,  | )**2-w(i,  J )**2) 
t = ak9**2  - ok  1 0 

if  (t.ge.O.)  go  to  9-* 
negr  = negr  1 

el  ( i » J ,2)  = 1.01*el( i*J , l)*sqrt (aklO > 
go  to  93 

uli,J)  = ( ak9+ sqr T ( t ) ) * 0. 5/akll 
r ( i » ] ) = el(i,J  ,i)/u(i,J  > 

p(itj)  = r(i,)l4(l.-u(i,))M2-vli.|l**2-w(i.|)**2) 

t = p (i , J ) 

if  (t.gt.O.)  go  to  36 

negp  = n e g p ♦ 1 

cont  inue 

return 

end 


subroutine  decode(l) 

dimension  p ( 42  * 42 ) *r ( 42  * 42 ) ,u  (42,42  )*v(42*42)  * w ( 42 * 42 ) 
dimension  pi  (42 , 42 ) * r I (42*42) * u I (42*42) * v I ( 42  » 4 2 ) * w I (42*42) 
dimension  e(42*42*4)  * f (42*42*4)  » g ( 42  * 4.2^,  4 ) * h (42*42*  4) 
dimension  el(42»42*4)  : 

common  /cl/  p»r*u*\/*w 

common  /c 7/  i i i « J J | * i i i I ♦ ] ) ) I * ) nl» ) t * 1 1 1 

common  /ell/  ak 

common  /cl4/  el,e,f * y*h 

common  /cl7/  negr, negp 

common  /cl  6/  pi  * r I * u I * v I * w I 

akll  - 1 . - ak 

common  A.2  0/  J a*  ) b 

do  36  i = l»iii 

do  36  ) = Ja.  I 1 J 

if  ( I . gt  .0  ) go  to  92 

vl (i* J ) = e( i , ] , 3) /e ( i * J *1) 

wf ( i , ] ) = eC i, J ,4)/e ( i , J ,1) 

ak9  = e ( i , J , 2 ) /e ( i * ) , 1) 

aklO  = 4.*akll*akMl.-vl (i,J >**2-wl (i, ])**2) 
t = ak9**2  - aklO 

if  (t*ge.0«)  go  to  37 
negr  - negr  + 1 

e( i * ) *2)  = 1.0i*e(i, ] , 1 ) * sor t ( ak 10 ) 
go  to  38 

u((i*J)  = (ak9+sqr t ( t ) ) *0 . 5/akil 

r il  ( i , J > = e(i,  J , 1)  /ul  (i , J ) 

Pl(i*)>  = rl (i,|)*(l.-ul (i, 1 Ci,|)**2-nl (it ])**2) 

t = pi (i, J) 

if  (t.gt.O.)  go  t o36 
negp  = negp  «•  1 

go  to  36 

v (i  , J ) = e (i  , J ,3) /e(  i , | ,1) 
w(i  , ) ) = e (i  ,) *4>/e<  i , ] ,1) 
ak9  = e (i  ,) ,2)/e(  i *|  ,1) 

aklO  = 4.*akll*ak* (l.-v(i, ) )**2-w  (i *J >**2> 

t = ak9**2  - aklO 

if  (t.ge.O.)  go  to  94 

negr  = negr  + 1 


51 


:-'r'  7 rJ!m  ~ 


go  to  93 

94  u(i»J)  = < ak9*sqrt  (t I )*Q.5/akll 

r < i * J > = e(i»)*l)/u(i,J) 

p ( i , 1 ) = r (i,J  )*(l.-u(i,J  >**2-v/(i,J  l**2-w(i,J  >**2> 
t = P (i,l  > 

it  (t.gt.O.)  go  to  36 
negp  = negp  ♦ 1 

36  continue 

return 
end 
c 
c 
c 

subroutine  prd 

dimension  e ( 42 , 42 ,4)  , f (42,42,4)  «g(42,42,*)  ,h(42,42»4) 

dimension  el(42»42*4) 

common  /c6/  J J n( a , xi u , t bb , e tr 

common  /c7 / ii i » 1 ] ] , i i i 1 , 1 ] ] I , I nl, ] t , ) tl 

common  /c9/  Jn 

common  /cl4/  el«e,f,g,h 

common  /ci5/  dz,al,bet«z 

common  /c2C/  la,  lb 

do  98  k = 1,4 

do  27  1 = Ja,  J 1J 

do27i-2,iii  ; | 

el ( i , J , k ) = e ( i , I , k)  - a I * ( f ( i , | ,k>  - f ( i'«  1 -1 , k > ) -be  t * ( g ( i + 1 , j ,k> 

& -g(i»J»ki)  - (h ( i« J ,k ) +2 . *e ( i , 1 ,k ) ) *dz 

27  continue 

do  97  1 = la, It 

97  el ( 1 , 1 , k ) = e(l,J,k)  - a I * < f ( 1 , J ,k  > -f  < 1 , ] -1 . K > > -be  t * ( g < 2 , 1 , k > 

1 -g  (1 , J , k) ) - Ih  ( 1,  J ,k)+2.*e(l,l,k))*dz 

do  201  1 = 1 tl,  J)) 

201  el ( 1 , 1 , k ) = e( 1, I ,k) 

98  continue 
return 
end 

c 

subroutine  cor 

dimension  e(42,42,4) ,f 142,42,4) ,g(42,42,4) ,h(42,42, 4) 

dimension  el(42, 42,41 

common  /c6/  J J n, a, xi u , tbb , e tr 

common  /c7/  i i i , J 1 1 , i i i I , 1 | J 1 , 1 nl , J t , 1 1 1 

common  /cl4/  el,e,f*g»h 

common  /cl5/  dz,al,bet,z 

common  /c2  0/  J a, Jb 

do  200  k = 1,4 

do  28  1 = la,  m 

do  28  i=  2 , i i i 

e ( i , 1 , k ) = 0.5* <e (i , ] ,k) *el f i, 1 ,k)-al* ( f C i, J +l,k)-f < i, J ,k) ) 

1 -be  t*(g(i,],K)-g(i-l,J,k))-(h(i,J,k)42.*el(i,J,k))*dz) 

28  continue 

do  99  J = | a, J t 

e < 1 , 1 , k ) = 0. 5* <e(l , 1 ,k) +el< 1, J ,k)  -a  I* ( f (1, 1 *l,k) - f ( 1, ) ,k) ) 

1 -bet*  (g(  2,j  ,k>-g(l,),k>>-(h(l,),k>*2.*elll,l,k)>*dz> 

99  continue 

200  continue 

return  ^ 

end 

subroutine  eo( I ) 


52 


dimension  p( 42,42), r ( g2 » 4 2)  ,u(42,42<»v(42,42)  , h (42 » 42  ) 

dimension  p(  (42, 42) , r I ( 42,42)  ,u)  (42, *2)  , v!  (42 , 42)  , * I ( 42 , 42) 

dimension  e(42,42,4)  , f(42,42,4),g(‘+2,42«4),h(42,42«4) 

dimension  el(42,42,4) 

common  /cl/  p,r,u,v,w 

common  /c4/  im, J m, im I , J mj 

common  /ell/  aK 

common  /cl4/  el,e,f»g,h 

common  /c!8/  p I , r 1 , u I , v I , w I 

common  /c2  0/  J a,  J b 

do  90  i = 1, im 

do  90  1 = Jb, | m 

if  ( I • gt  .0  ) go  to  91 

e(i,),l)  = r I ( i , ] ) *u  I ( i , J ) i 

e(i,j,2)  = p I ( i , ) ) *ak  + e ( i , ] ,1 ) *u 1 ( i , J ) 
e(i,],3)  = e ( i , J , 1) * v I ( i , J ) 
e(i,J,4)  = e ( i , ] , 1) * w t ( i , J > 
go  to  90 

91  e ( i , ] , 1 ) = r ( i , J )*u(i,J  ) 

e ( i , ] , 2 ) = p ( i , ) ) * aK  + e { i , ) , i > *u( i , J ) 
e(i, J ,3)  = e (i , J , 1 )*v ( i, ] ) 
e(i»l,4)  = e ( i , 1 , 1) * w ( i , J ) 

90  continue 

return 

end  •:  ‘ '!  !‘  ' 

c 

subroutine  fghU) 

dimension  p(42,42),r(42,42) ,u(42 *42 )*v( 42,42)  ,w(42,42) 

dimension  pi  ( 42 , 42 ) « r I ( 42 , 42  ) , u 1 (42 , -*2  ) , v I (42 , 42 ) , w I (42,42) 

dimension  et(42*42),xi(42,««2)  , ch i ( 42 ) , psi  ( 42 ) 

dimension  et)(42,42),xil(42,42),chi  l(42),psil‘4  2) 

dimension  el‘,2,42,4)  , f (42,42,4)  , g ( 42 , 42 , 4)  ,h«  ,2,42,4) 

dimension  el(42,42,4) 

dimension  fb(4),gb(4) 

common  /cl/  p,r»u,v,w 

common  /c4/  im, ) m, im J , J ml 

common  /c6/  ) ) n, a,xi u, tbb,etr 

common  /c8/  x i , e t , xi I , e t I 

common  /clQ/  xi r , chi , ps  i ,ch i ( , ps i 1 

common  /ell/  aK 

common  /cl  3/  s 

common  /cl4/  el,e,f»g,h 

common  /cl6/  de I chib , de I ps i b , de 1 ch i I 

common  /cl  8/  p I , r I ,u I ,v  I , m I 

common  /c 2 0/  J a,  J b 


do  34  i 
do  34  J 
if  ( I .gt 
v v 

MM 

f b ( 1 ) 
f b ( 2 ) 
f b ( 3 ) 
f b ( 4 > 
gb  ( 1 ) 
gb  ( 2 ) 
gb  ( 3 ) 
gb  ( 4 ) 
do  05 


K - 


f (i t J ,k> 


et I (i, ) ) *ul ( i , ) > 
xi  I ( i, J >*u  I ( i , J ) 


= 1,  im 
= 1 b,  1 m 
.0 ) go  t o 86 
= vl (i , J )/s  - 
= w I ( i , J ) / s - 
= vv*r I ( i , ] ) 

= vv*r I ( i , J ) * u I ( i , J ) 

= vv*r I ( i , J ) * v I ( i , ) ) 

= vv*r  I (i,  ))*m)  li,J  ) 

= ww  *r  I ( i , J ) 

= mm*t I ( i , ) ) *u I ( i , | ) 

= MW*r  I ( i , J ) * v I ( i , ) ) 

= ww*r l(i,])*wl(i,]) 

1,4 

= (chi  I ( J ) *gb(k) tf b (K>  >*chi I ( J )/etl (i,J) 


ak*pl ( i , ] ) *et(  (i,J) 
ak*pl ( i , ] )/s 


aK*p  I ( i , | ) 4xi I ( i , ) ) 
aK*  p I ( i , J ) /s 


53 


SM  i , J , h ) 

, - ! \ 1 - i i » ) r i i.>  « • , i 

h( i , ) , h> 

ltooMb(K)«-(l*+2.;cnil  l J )¥lDo)*gu(n))/(psil  (i)-xiu) 

f ( i • 1 « K ) 

= 

f < i « | , k) /de ichi 1 

85 

g ( i , 1 , K) 
go  to  34 

= 

g ( i , J , k) /psi 1 ( i m 1 ) 

86 

vv 

= 

v(l,J)/s  - et ( i , ) ) *u ( i f J ) 

MW 

= 

w(i,J)/s  - xi ( i , ] ) *u  (i « |) 

f b<  1) 

= 

vv*r  ( i , ] ) 

f b(2> 

= 

vv*r  ( i t J ) *u ( i » ) ) - ak*p  <i  , J ) *e  t ( i , J ) 

f b ( 3 ) 

vv*r  ( i , J ) * v ( i * ] ) * ak*p(i,J)/s 

f b ( 4 ) 

v»/*r  ( i » ) ) *w  ( i , J ) 

gb  ( 1 ) 

= 

ww¥r  ( i » J ) 

gb  ( 2 > 

= 

ww*r  ( i , ) ) *u  ( i , ) ) - ak*p ( i , ) ) *xi ( i , J ) 

gb(  3 ) 

= 

ww*r  ( i , ) ) ( i , | ) , 

gb  ( 4 ) 

= 

ww’rlii] )¥w(ii))  + ak*p(i,J)/s 

do  35  k = 

If 

4 

f (i,),k> 

= 

(chi(J)*gb(k)  * fb(k))*chi(j)/et(i«j) 

g(i  * j »k) 

= 

(psi ( i ) * f o (k)  ♦ gb ( k ) ) * ps i ( i ) / ( xi  ( i , | ) - xi r ) 

h(i , l ,k) 

(psi(i)*f(i,]«k)+chi(])*g(i»J,K))/(l.-osi(i)*chi(j)> 

&> 

-tb(k)*psi (i)/(xi (i,J ) -xi r) 

8. 

-gb(k)*chi(J)/et(i,i  ) 

f ( i » ) * k ) 

=r 

f ( i » J * k) /de 1 chi b 

g(i » J * k) 

g(  i , J « k) /de 1 psib 

35 

cont inue 

34 

cont inue 

return 

end 

; , * 1 1 ! • ' 

c 

c 

c 

subroutine  pr 

dimension  p(42»42),r(42,42),u(42,42),v(42,42)  , m ( 42 , 42 ) 
dimension  e t ( 42 ♦ 42 ) » x i ( 42 , 42 ) 

dimension  pi  (42,42), rl  (42,42)  , u I (42 , 4 2 ) , v I ( 42 , 42 ) « w I (42 ,42) 

dimension  e 1 1 ( 42 , 42)  , xi I ( 42 , 42) 

common  /cl/  p,r,u,v,w 

common  /c4/  i m , J m , im  I , J m I 

common  /c7/  i i i * ] | J * i i i I , J ) J I » 1 nl» J t , J 1 1 

common  /c9/  Jn 

common  /clfi/  pi »r I »ul  ,vl * w I 


c partial  output 


print 

110 

print 

113* 

<p(l, J ),J-l,)tl,2> 

print 

114, 

(pi (1, J 1) , J 1 = 1, J 1 1 , 2 ) 

print 

148, 

(rl (1, J 1) , ) 1 = 1, J tl,2) 

print 

149, 

(vt  (1,)  1) , ) 1 = 1, J 1 1 , 2 ) 

print 

11C 

print 

113, 

( p ( i , 2 ) , i=l , im, 2 ) 

print 

114, 

(pi (i 1,2) , i 1 = 1, im  1,2) 

print 

110 

do  40 

k = 1 

,3 

40 

print 

113, 

( p ( k, | ) ,J=Jtl,Jm) 

print 

110 

do  63 

k = 1 

,3 

63 

print 

114, 

( d 1 (k, J 1)  , ) l = J ti, Jml ) 

pri  nt 

110 

return 

c 

110 

format 

( 2x ) 

| 

113 

format 

C‘P 

" ,2x,20f8.4) 

54 


14  8 

149 

c 

c 

c 


c 


c 

ill 

121 

c 

c 

c 


c 

52 

53 

46 


r\  . . . ,2x>~Gf8«4) 

I or  mat  C‘r  I" , 2x , 20  f 8 . 4 ) 

format  ("v  I",  2x,  2 0 f 8 « 4) 

end 


subroutine  prcoord 
dimension  e t (42 , 4 2) , x i ( 42 , 42 ) 
dimension  et  1 ( 42 , 42)  , x i I ( 42 , 42) 
common  /c4/  im , | m , im I , ) m I 
common  /c3/  xi , e t , xi I , e t I 
common  /c9/  )n 

partial  coord  output 
J nl  = J n 4-  i 
print  121 

print  111,  < et (1, | ) , 1 =1, | nl, 4) 

print  121 

return 

format  (**e  t*' , 2x , 20  f 8 • 4 ) 
format  (Ex) 

end 


subroutine  prf in 

dimension  p(42,42),r(42,42),u(42*42),v(42*42)  , w ( 42,  42 ) 

dimension  et(42,42l,xi(42»42),chi(42)  , ps  i ( 42 ) 

dimension  pi  ( 42 , 42 ) , r I (42,42)  ,ul  (42,42)  , v I (-42,42)  , w I (42,42) 

dimension  etl ( h2 , 42) , xi I ( 42 , 42 ) , chi  1 ( 42 ) , ps  i I (42) 

dimension  etb ( 42 ) , xi b ( 42 ) 

dimension  p2 ( 42 , 42) , p3 ( 42, 42 ) 

common  /cl/  p,r,u,v,a 

comnon  /c2/  ga , a m , be , pi , amd 

common  /c3/  hh , dec , d cp , tde , occ J , dep I 

common  /c4/  im, ) m, im I , J ml 

common  /c6/  J J n , a , xi u , t bb , e t r 

common  /cl / i i i , ] J J , i i i I , ] ) ] I , ) n 1, J t , ] 1 1 

common  /ce/  xi , e t , xi  l , e 1 1 

common  /ccJ/  Jn 

common  /c 10/  xir , chi , psi ,chi I , ps i 1 
common  /cl8/  p I , r I ,u I , v I , w I 
coordinates 

print  125 
do  52  i=l» im,2 

print  126,  ( et ( i , J ) , J =1 , J m, 2 ) 
print  125 
do  53  i=l, im»2 

print  127,  ( xi  ( i , J ) , J =1 , J m*  2 ) 

do  46  J = J n , J m 

etb(J)  = (chi(J)*xiu)/(l.+cni())*tbb) 
x i b ( J ) = tbb4  e tb ( ] ) 

print  125 

print  126,  ( etb  ( 1 ) , 1 ■=  In, ) m) 


print  129,  ( xi  o ( J ) , J = J n , J m) 

print  125 

do  64  i I = 1,  iml , 2 


64  print  138,  ( et I ( i I , J I ) , J I = 1 , J m I , 2) 

print  125 


d o 65 

i 1 = 

1,1ml  ,2 

65 

print 

139, 

( xi  l ( i i , ) 1 ) , J 1 = 1 , 1 ml , 2 

print 

125 

c 

field 

ak9  - 

: 2 ./ (ga* (am**2 1 ) 

do  47 

] = :. 

J m 

do  4 7 

1=1, 

i m 

47 

p2 ( i , ] 

) = 

p(i,J )/p(im,Jm) 

do  9 5 

i=  1 

, i m 

95 

pri  nt 

130, 

(p2(i , | ) , 1 =1,14) 

print 

125 

do  2 02 

: i = 

1 , im 

202 

print 

130, 

(p2 (i , ] ) , ] =15,28) 

print 

125 

do  5 4 

i=  1, 

54 

print 

130, 

(p2 ( i , J > , ) =29, ] m) 

pri  nt 

125 

do  5 5 

i- 1» 

im,6 

55 

pri  n t 

131, 

( r ( i , J ) , ] =1 , J m,4 > 

print 

125 

do  56 

i = 1, 

im,6 

56 

pr  i n t 

132, 

( u( i , ] ) , J =1, Jm,4) 

print 

125 

do  57 

i=  1, 

1 m ,6 

57 

print 

133, 

( v (i , J ) , 1 =1, J m,4) 

print 

125 

do  53 

1=1, 

im,6 

58 

print 

134, 

( w( i , ] ) , 1 =1, | m,4) 

pri  n t 

125 

J ln4 

= 1 n - 4 

J Jn2 

= J n + 2 

pri  n t 

131, 

( r ( 1 , J ) , j = J J n4,  J J ri2 > 

print 

132, 

( u (1 , J ) , ) = J ] n4, ) ] n2) 

pri  n t 

133, 

( m (1, ] ) , J =] 1 n4 , 1 | n2) 

print 

134, 

1 w(l, J > , J = J J n4, ]] n2) 

print 

125 

do  66 

J 1 = 

1 , ] m 1 

do  66 

i 1 = 

1 , i m 1 

66 

p3  I i 1 , 

1 l> 

= pi  (i 1 , J ») /pi ( i m 1 , ) ml  ) 

do  96 

i 1 = 

1 , i m 1 

96 

print 

140, 

(p3< i 1 , 1 II , j 1=1,14) 

print 

125 

do  2 0 3 

i 1 

= 1 , i m 1 

203 

print 

14  0, 

(p3 (i 1 , J 1 ) , | 1=15,28) 

print 

125 

do  67 

i 1 = 

1,2 

67 

pr  i n t 

140, 

(p3 ( i 1 , ] 1 ) , ) 1=29, ] ml ) 

print 

125 

do  6 8 

i 1 = 

1 , i m 1 ,6 

68 

print 

141, 

(rl  (i 1 , J 1 ) , J 1 = 1, J ml, 4) 

print 

125 

do  69 

i 1 = 

1 , i m 1 ,6 

69 

print 

142, 

(ul 111 , | 1 > , J 1 = 1, J ml ,4) 

print 

125 

do  7 0 

i 1 = 

1,  iml  ,6 

70 

print 

143, 

(vl (i 1 , | 1 ) , 1 1 = 1,1  ml, 4) 

print 

125 

do  71 

i 1 = 

1 , i m 1 ,6 

71 

print 

144, 

(wl  (i  1 , J 1 ) , ] 1 = 1, J ml ,4) 

print 

125 

1 


56 


c I p =0. 
c I p I = C. 

J Jn3  = It  - 1 
do  60  1 = 2,  J ] ni 

60  clip  = cl P ♦ (p2 (1, ) *1) +p2(l , | ) )* (et (1, )*1) -el  (i, J ) ) /2. 

clp  = clp  ♦ (p2  ( 1,2)  +p2(l,l ) )*  (et  ( i.,  2 J-et  (1, 1 ) >/4. 
ctp  = clp  + p2(l,JJn)Mi.-et(l,)Jn)> 
do  7 4 J I = 2,  J J n3 

74  clpl  = clpl  + (p3(l,  J I +1)  +p3  (1,  ] I ) > Met  I (1,  ] I +l)-et  I (1,  J I ) ) /2. 

cDpI  = clpl  ♦ (p3 (1,2) +p3(l,l ) )* (et I (1,2) -etl ( l, 1) )/4. 
cdpl  = clpl  + p3  ( 1,  ] | n) *(l.-et I ( l, ) ] n) ) 
c 8 = ak9  v (c I p-c  J pi  I 

cd  = ak9*clp*Tde 
c I d = c 1/  cd 
print  137,  cl,cd,cld 
print  125 
do  48  J = 1 , ) m 
do  48  i - lfim 

48  p2(i,J)  = ak9* ( p2 ( i * ) ) - 1 . ) 

print  136,  ( p2 ( 1 * 1 ) , 1=1,14) 

print  13  b,  (p2 ( 1 , J ) , ) =15,2 3) 
print  136,  ( p2 (1 , | ) , ) =29, J m) 

print  125 

do  72  J I = 1,  | ml  •;  ' " ' 

do  72  i I = 1 , i ml 

72  p3  ( i I , ] I ) = ak9M  p3(  i I , | I ) -1 . ) 

print  1 1, 5,  (p3  (1 , ) ) , J =1 , 12 ) 

print  145,  (p3 (1 , J ) , ] =13, 1 m) 

print  125 
return 
c 

125  format  (2x) 

126  format  ( "e t" , 2x, 16 f 8 . 4 ) 

127  format  ( *'x  i”,  2x,  16  f 3 . 4) 

128  format  ( "e tP“ , x, 16 f 9 . 4) 

129  format  ( *'x  id** , x,  16  f 8 . 4) 

130  format  ( “p*' , 3x  , 1 6 f 8 . *♦ ) 

131  format  ( “r  **,  3x,  16  f 3.  4 ) 

132  format  ( "u", 3x, 1 6 f 6. 4 ) 

133  format  ( ‘V,  3x,  1 5f  8.  4) 

134  format  ( 'V' , 3x , 16 f 8. 4 ) 

136  format  (*'cp,',2x,16f8.4) 

137  format  ( **c  I =** , f 8 . 4,  hx,  "cd  =" , f 8 . 4 , 4 x ,”c  I /c  d =**,f8.4) 

138  format  ( "e t 1 M, x, 16 f 8 . 4) 

139  format  ( *'x  i I'*,  x , 16  f 8 . 4 ) 

140  format  ( *'p  I"  , 2x,  16  f 8 . 4 ) 

141  format  ( V I " ,2x, 16 f 3 . 4 ) 

142  format  ( "u I “ , 2 x, 16 f 8 . 4) 

143  format  ( "v  I *' , 2x , 16  f 8 . 4 ) 

144  format  (*'w  I"  ,2x,  16f  8 . 4) 

145  format  ( *'c  p I ’*,  x,  16  f 8 . 4 ) 

c 

end 

c 

c 

subrout  i re  mod  ( I ) i 

dimension  p(42,42),r(42,42),u(42,42),v(42,42) ,w(42»42) 

dimension  pi  ( 42 , 42 ) , r I (42,42)  ,ul  (>42,42)  *vl  (42,42)  ,rtl  (42,42) 


^ — -■'-■'T 


82 

c 

83 


81 

c 


d i mens 

common 

common 

common 

common 

common 

common 

common 

common 

common 

common 

if  (l 

do  82 

psi  s 

di 

1 i 

d xi  r 
Pfl*)  I 
r « 1 , J I 
utl,  J I 
V ( 1 » J I 
W ( 1 , J i 
return 


57 

i cn  e t l ( 2 , h 2 ) < x i I ( 4 l < 4 L t * ch  i I ( 4 2 ) > .>  s i 1 ( 4 c ) 

/cl  / p i r , u i v,  w 

/c 3/  hh , dec ,d cp « tde» dec  I * dep  I 

/c4/  iin>]m*imt,]ml 

/c6/  ] ) n » a , xi u , t Db » e tr 

/cl / i i i i J | ] t i i i I , J 1 ] I . ] nl, ) t , ] 1 1 

/c8/  xi,et,xil,etl 

/c 9/  ] n 

/clQ/  x ir t chi , ps i t ch i I , ns i I 
/cl  8/  pi ,rl ,ul  tvl i«l 
/cl  3/  s 

g t .0  ) go  to  83 

] I = J tl,  J ] J 

- tde*(l.+(tbb-xiu)*chil(J  l))/<sMl.+hh*chil{]l))) 
= 1.  +■  (psi  s/dcp  I ) 

= if  i x ( d i ) 

(xil  (i  I , 1 I ) - x i ( 1 * J I ) ) / ( x i I ( i I , J I ) -x  i I ( i I +1 t J I )) 

= pl(il,]i)Mpl(il+lf)l)-pl(il»)l)>*dxir 
= rl(il»JI)+(rl(il+lfJI)-rl(il,]l))^dxir 
= ul(ilt]l)MuHi]  + l,]l)-ul(iltJ3>)*d*ir 
= vl(il,JI)Mvl{i1  + lt|l)-vl(il,]ll))»dxir 
) = wl(il,jl)+(wl(il+l»|l)-wl(il,)l))*dxir 


do  81  1 
ps  i s = 
di 

i = 

dxi  r = 
Pi  (It  ) ) 
rl  ( 1,  J ) 
uJ  ( 1,  | > 
vt  (1,  J > 
wC (1,  J ) 
return 

end 


= JfltlJ  1 

((xiu-xir)*tOD*chi(J  )-xir)  / ((etr*tpp-xiu)*chi ()  l+etr) 
1.  + ( ps  i s+h  h ) / dc  > 
if  i x (di  ) 

(xil  (ltJ  > — x i ( i * ) ))/(xi(i+ltJ ) -xi ( i » J ) ) 

= p ( i,  ) ) + (o(  i+lt 1 )-p ( i t ) )> * cxir 
= r(i,l>  + (r(i  + lt|)-r(itj))*axir 
= u(itl)  + lu(iti,))-uli»)))^cxir 
= v(  i t ) ) + (v/(  i +1 , ) ) -v  ( i t J ) ) * axir 
- w(i»])+(w(i+ltj)-w(itj))*dxir 


Normal  Section  Angle  Dependence  on  Span  and  Shock 
Angle  for  Design  Mach  Numbers  2 and  4 


Table  2 


On-Design  Normal  Section  Shock  Angle 


Design  Attached  Shock  Conditions 
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Table  4 

Case  Constants 


Attached  Algorithm 


s* 

M 

00 

x jm 
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Jm 
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Detached  Algorithm 


s* 

M 

00 

^m 

X 

In 

nR 

nM 

r 

'n 

^C 

1 

1.6 

24 

x 

24 

1. 

35 

1. 

30 

1. 

05 

2.637 

2.276 

2 

1.7 

28 

X 

24 

1.467 

1.' 

40 

5.882 

1.225 

1.5 

24 

X 

24 

0.980 

1.4 

24 

X 

24 

1.225 

1.35 

16 

X 

24 

0.980 

1.3 

20 

X 

24 

1 

1 

0.980 

Figure  12.  PRESSURE  DISTRIBUTIONS  FROM  SURFACE  TO  FREE  STREAM  ALONG 
SYMMETRY  PLANE:  s*=2,  BD=35° 

(a)  6=  5.75°,  M = 1.6,  1.8,  2.0 

* oo 


0.24 


Figure  13.  SURFACE  PRESSURE  COEFFICIENT  DISTRIBUTION;  s*=2,  6n=35 


1.0  1.2  1.4  1.6  1.8  2.0 


Figure  16.  DRAG  COEFFICIENT  VARIATION  WITH  MACH  NUMBER; 
6 = 5.75°,  Bd  = 35° 
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